#include "allcode.h"

/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja                                  *
 *   harish.kukreja@gmail.com                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//
// C++ Implementation: physiksvector
//
// Description: Definition of PhysiksBase::Vector.
//
//
// Author: Harish Kukreja <harish.kukreja@gmail.com>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//



PhysiksBase::Vector::Vector(unsigned int size, double init) : dimension(size) {
    cont = gsl_vector_alloc(size);
    gsl_vector_set_all(cont, init);
}

PhysiksBase::Vector::Vector(const gsl_vector *init) : dimension(init->size) {
    cont = gsl_vector_alloc(init->size);
    gsl_vector_memcpy(cont, init);
}

PhysiksBase::Vector::Vector(const Vector &init) : dimension(init.dimension) {
    cont = gsl_vector_alloc(init.dimension);
    gsl_vector_memcpy(cont, init.cont);
}

PhysiksBase::Vector::~Vector() {
    gsl_vector_free(cont);
}

void PhysiksBase::Vector::scale(double fact) {
    gsl_vector_scale(cont, fact);
}

void PhysiksBase::Vector::set_zero() {
    gsl_vector_set_zero(cont);
}

void PhysiksBase::Vector::set(unsigned int index, double value) {
    gsl_vector_set(cont, index, value);
}

void PhysiksBase::Vector::set(const Vector &copy) {
    gsl_vector_memcpy(cont, copy.cont);
}

double PhysiksBase::Vector::norm() const {
    return gsl_blas_dnrm2(cont);
}

PhysiksBase::Vector PhysiksBase::Vector::normalize() const {
    Vector res(*this);
    gsl_vector_scale(res.cont, 1.0 / gsl_blas_dnrm2(res.cont));
    return res;
}

PhysiksBase::Vector PhysiksBase::Vector::duplicate() const {
    return Vector(*this);
}

unsigned int PhysiksBase::Vector::dim() const {
    return dimension;
}

double PhysiksBase::Vector::get(unsigned int index) const {
    if (index < cont->size && index >= 0) return gsl_vector_get(cont, index);
    else return 0;
}

PhysiksBase::Vector PhysiksBase::Vector::cross(const PhysiksBase::Vector &rhs) const {
    PhysiksBase::Vector res(3, 0);
    res.set(0, this->get(1) * rhs[2] - this->get(2) * rhs[1]);
    res.set(1, this->get(2) * rhs[0] - this->get(0) * rhs[2]);
    res.set(2, this->get(0) * rhs[1] - this->get(1) * rhs[0]);
    return res;
}

PhysiksBase::Vector PhysiksBase::operator+(const PhysiksBase::Vector &lhs, const PhysiksBase::Vector &rhs) {
    PhysiksBase::Vector res(lhs);
    gsl_vector_add(res.cont, rhs.cont);
    return res;
}

void PhysiksBase::operator+=(PhysiksBase::Vector &lhs, const PhysiksBase::Vector &rhs) {
    gsl_vector_add(lhs.cont, rhs.cont);
}

void PhysiksBase::operator-=(Vector &lhs, const Vector &rhs) {
    gsl_vector_sub(lhs.cont, rhs.cont);
}

PhysiksBase::Vector PhysiksBase::operator-(const Vector &lhs, const Vector &rhs) {
    Vector res(lhs);
    gsl_vector_sub(res.cont, rhs.cont);
    return res;
}

PhysiksBase::Vector PhysiksBase::operator-(const Vector &rhs) {
    Vector res(rhs);
    gsl_vector_scale(res.cont, -1.0);
    return res;
}

double PhysiksBase::operator*(const Vector &lhs, const Vector &rhs) {
    double res;
    gsl_blas_ddot(lhs.cont, rhs.cont, &res);
    return res;
}

PhysiksBase::Vector PhysiksBase::operator*(const Vector &lhs, double fac) {
    Vector res(lhs);
    gsl_vector_scale(res.cont, fac);
    return res;
}

PhysiksBase::Vector PhysiksBase::operator*(double fac, const PhysiksBase::Vector &rhs) {
    PhysiksBase::Vector res(rhs);
    gsl_vector_scale(res.cont, fac);
    return res;
}

double PhysiksBase::Vector::operator[](unsigned int k) const {
    return gsl_vector_get(cont, k);
}

gsl_vector *PhysiksBase::Vector::core() const {
    return cont;
}


/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja   *
 *   harish.kukreja@gmail.com   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//
// C++ Implementation: collision
//
// Description: Definitions of collision detection & reaction functions.
//
//
// Author: Harish Kukreja <harish.kukreja@gmail.com>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//


using PhysiksBase::Computation::compute_distance;
using PhysiksBase::Computation::compute_connecting_vector;

// Collision detection functions //////////////////////////////////////////////
bool PhysiksBase::Collisions::detect_collisions(const World &world) {
    for(uint i = 0; i < world.num_objects(); i++) {
        // Perhaps one could do a pairwise comparison, by calling one of the detect_* functions
        // depending upon each object's type.
    }
    return false;
}

///////////////////////////////////////////////////////////////////////////////
bool PhysiksBase::Collisions::detect_point_sphere(const PhysicalObject &obj, const Sphere &sph) {
    return compute_distance(obj, sph) < sph.get_radius();
}

///////////////////////////////////////////////////////////////////////////////
bool PhysiksBase::Collisions::detect_point_plane(const PhysicalObject &obj, const Plane &plane) {
    /** Find the minimum distance between the point and the plane. If the `base' of the plane-point
     * interconnect lies in the bounded region, return true. */
    Vector temp(obj.position - plane.position);
    Vector distance(plane.normal() * temp * plane.normal());
    if (!plane.bounded()) {
        return distance.norm() < 1.0E-3;    // This arbitrary value should not be hard-coded.
    } else {
        if (plane.contains(obj.position - distance)) return distance.norm() < 1.0E-3;
        else return false;
        // This doesn't take into account "skimming" collisions, where only a small portion of the
        // sphere may strike the edge of the bounded plane.
    }
}

///////////////////////////////////////////////////////////////////////////////
bool PhysiksBase::Collisions::detect_sphere_sphere(const Sphere &sph1, const Sphere &sph2) {
    return compute_distance(sph1, sph2) < (sph1.get_radius() + sph2.get_radius());
}

///////////////////////////////////////////////////////////////////////////////
bool PhysiksBase::Collisions::detect_sphere_plane(const Sphere &sph, const Plane &plane) {
    Vector temp = sph.position - plane.position;
    Vector distance = plane.normal() * temp * plane.normal();
    if (plane.bounded()) {
        // How does one account for bounded planes? Hmmm....
        // This doesn't take into account "skimming" collisions, where only a small portion of the
        // sphere may strike the edge of the bounded plane.
        if (plane.contains(sph.position - distance)) return distance.norm() < 1.0E-3;
        else return false;
    } else {
        /** Find the minimum distance between the centre of the sphere and the plane, and compare
         * it to the sphere's radius. */
        return distance.norm() < 1.0E-3;    // This arbitrary value should not be hard-coded.
    }
}

//=============================================================================

// Collision reaction functions ///////////////////////////////////////////////
void PhysiksBase::Collisions::react_point_sphere(PhysicalObject &obj, Sphere &sph) {
    obj.net_force += compute_connecting_vector(obj, sph) * sph.get_elasticity();
}

void PhysiksBase::Collisions::react_point_plane(PhysicalObject &obj, Plane &plane) {
}

void PhysiksBase::Collisions::react_sphere_sphere(Sphere &sph1, Sphere &sph2) {
    double dist = compute_distance(sph1, sph2);
    Vector temp = (sph2.position - sph1.position) * (1 / dist);
    sph1.net_force += temp * sph2.get_elasticity() * (sph2.get_radius() - dist);
    sph2.net_force -= temp * sph1.get_elasticity() * (sph1.get_radius() - dist);
}

void PhysiksBase::Collisions::react_sphere_plane(Sphere &sphere, Plane &plane) {
}
/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja                                  *
 *   harish.kukreja@gmail.com                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//
// C++ Interface: computations
//
// Description: Definitions of miscellaneous computation functions used in Physiks. Also includes
//              definitions for functions that modify the states of PhysicalObjects according to
//              fields.
//
// Author: Harish Kukreja <harish.kukreja@gmail.com>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//



using namespace std;

using PhysiksBase::Vector;

void PhysiksBase::Computation::apply_force_fields(vector<PhysicalObject *> &objTable,
        const vector<ForceField *> &force_fields) {
    for (uint i = 0; i < force_fields.size(); i++)
        for (uint k = 0; k < objTable.size(); k++)
            objTable[k]->net_force += (*force_fields[i])(objTable[k]->position);
}

void PhysiksBase::Computation::apply_acceleration_fields(vector<PhysicalObject *> &objTable,
        const vector<AccelerationField *> &accel_fields) {
    for (uint i = 0; i < accel_fields.size(); --i)
        for (uint k = 0; k < objTable.size(); --k)
            objTable[k]->net_force += objTable[k]->mass
                    * (*accel_fields[i])(objTable[k]->position);
}

void PhysiksBase::Computation::apply_velocity_fields(vector<PhysicalObject *> &objTable,
                                                     const vector<VelocityField *> &vel_fields) {
}

double PhysiksBase::Computation::compute_distance(const PhysicalObject &a, const PhysicalObject &b)
{
    return Vector(a.position - b.position).norm();
}

Vector PhysiksBase::Computation::compute_connecting_vector(const PhysicalObject &a,
                                                           const PhysicalObject &b) {
    return Vector(b.position - a.position).normalize();
}

void PhysiksBase::Computation::update_forces(vector<PhysicalObject *> &objTable) {
    for (uint i = 0; i < objTable.size(); i++)
        for (uint j = i + 1; j < objTable.size(); j++) {
        objTable[i]->net_force += objTable[i]->mass
                * objTable[j]->gravitational_field(objTable[i]->position);
        objTable[i]->net_force += objTable[i]->charge
                * objTable[j]->coulomb_acceleration_field(objTable[i]->position);
        objTable[j]->net_force += objTable[j]->mass
                * objTable[i]->gravitational_field(objTable[j]->position);
        objTable[j]->net_force += objTable[j]->charge
                * objTable[i]->coulomb_acceleration_field(objTable[j]->position);
        }
}

double PhysiksBase::Computation::compute_gravitational_force(const PhysicalObject &a, 
                                                             const PhysicalObject &b) {
    double dist = compute_distance(a, b);
    return (GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT * a.mass * b.mass / (dist*dist));
}

double PhysiksBase::Computation::compute_coulomb_force(const PhysicalObject &a,
        const PhysicalObject &b) {
    double dist = compute_distance(a, b);
    return (a.charge * b.charge / (GSL_CONST_MKSA_VACUUM_PERMITTIVITY * 4.0 * M_PI * (dist*dist)));
}

void PhysiksBase::Computation::update_velocities(vector<PhysicalObject *> &objTable,
        double time_step) {
    for (uint i = 0; i < objTable.size(); i++)
        objTable[i]->velocity += (objTable[i]->net_force * (time_step / objTable[i]->mass));
}

void PhysiksBase::Computation::update_positions(vector<PhysicalObject *> &objTable,
        double time_step) {
    for (uint i = 0; i < objTable.size(); i++) {
        objTable[i]->position += (objTable[i]->velocity * time_step);
        objTable[i]->position += (objTable[i]->net_force * (time_step * time_step /
                (2 * objTable[i]->mass)));
    }
}

void PhysiksBase::Computation::set_force_vectors_zero(vector<PhysicalObject *> &objTable) {
    for (uint i = 0; i < objTable.size(); i++) objTable[i]->net_force.set_zero();
}
/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja                                  *
 *   harish.kukreja@gmail.com                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//
// C++ Implementation: fields
//
// Description: Implementations of in-build fields.
//
//
// Author: Harish Kukreja <harish.kukreja@gmail.com>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//


using PhysiksBase::PhysicalObject;
using PhysiksBase::Vector;

///////////////////////////////////////////////////////////////////////////////
PhysiksBase::NewtonianGravitationalField::NewtonianGravitationalField(PhysicalObject &obj) {
    center = &obj;
}

Vector PhysiksBase::NewtonianGravitationalField::operator()(const Vector &pos) const {
    Vector dist(pos - center->position);
    return Vector(GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT * center->mass / (dist * dist) * dist.normalize());
}

PhysiksBase::NewtonianGravitationalField::~NewtonianGravitationalField() {}


///////////////////////////////////////////////////////////////////////////////
PhysiksBase::CoulombAccelerationField::CoulombAccelerationField(PhysicalObject &obj) {
    center = &obj;
}

Vector PhysiksBase::CoulombAccelerationField::operator()(const Vector &pos) const {
    Vector dist(pos - center->position);
    return Vector(center->charge / (GSL_CONST_MKSA_VACUUM_PERMITTIVITY * 4.0 * M_PI * (dist * dist)) * dist.normalize());
}

PhysiksBase::CoulombAccelerationField::~CoulombAccelerationField() {}
/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja                                  *
 *   harish.kukreja@gmail.com                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//
// C++ Implementation: objects
//
// Description: Definitions of various physical Objects that are found in the Physiks interactions
//              environment.
//
// Author: Harish Kukreja <harish.kukreja@gmail.com>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//


using PhysiksBase::Vector;

// Definitions for PhysiksBase::PhysicalObject ////////////////////////////////

PhysiksBase::PhysicalObject::PhysicalObject() : position(3, 0.0), velocity(3, 0.0), net_force(3, 0.0), gravitational_field(*this), coulomb_acceleration_field(*this) {
    mass = charge = 0.0;
/*    position = new Vector(3, 0.0);
    velocity = new Vector(3, 0.0);
    net_force = new Vector(3, 0.0);
    
    gravitational_field = new NewtonianGravitationalField(*this);
    coulomb_force_field = new CoulombForceField(*this);*/
}

PhysiksBase::PhysicalObject::PhysicalObject(double m, double c, const Vector &pos, 
                                            const Vector &vel) : position(pos), velocity(vel), net_force(3, 0.0), gravitational_field(*this), coulomb_acceleration_field(*this) {
    mass = m;
    charge = c;
/*    position = new Vector(pos);
    velocity = new Vector(vel);
    net_force = new Vector(3, 0.0);
    gravitational_field = new NewtonianGravitationalField(*this);
    coulomb_force_field = new CoulombForceField(*this);*/
}

PhysiksBase::PhysicalObject::~PhysicalObject() {}

void PhysiksBase::PhysicalObject::draw() const {
/*    glBegin(GL_POLYGON);
        glVertex3f(position[0], position[1], position[2]);
    glEnd();*/
    glTranslatef(position[0], position[1], position[2]);
    glutSolidSphere(1.0, 180, 180);
}


// Definitions for PhysiksBase::ExtendedObject ////////////////////////////////

PhysiksBase::ExtendedObject::ExtendedObject() : PhysicalObject(), angular_velocity(3, 0.0) {
}

PhysiksBase::ExtendedObject::ExtendedObject(double mass, double charge, const Vector &position, const Vector &velocity) : PhysicalObject(mass, charge, position, velocity), angular_velocity(3, 0.0) {
}

PhysiksBase::ExtendedObject::~ExtendedObject() {}


// Definitions for PhysiksBase::Sphere ////////////////////////////////////////

PhysiksBase::Sphere::Sphere() : ExtendedObject() {
    radius = 1.0;
}

PhysiksBase::Sphere::Sphere(double mass, double charge, const Vector &position,
                            const Vector &velocity, double radius) : ExtendedObject(0.0, 0.0, position, velocity) {
    this->radius = radius;
}

PhysiksBase::Sphere::~Sphere() {}

void PhysiksBase::Sphere::draw() const {
    glTranslatef(position[0], position[1], position[2]);
    glutSolidSphere(radius, 180, 180);
}

double PhysiksBase::Sphere::get_elasticity() const {
    return elasticity;
}

void PhysiksBase::Sphere::set_elasticity(double val) {
    elasticity = val;
}

double PhysiksBase::Sphere::get_radius() const {
    return radius;
}

void PhysiksBase::Sphere::set_radius(double radius) {
    this->radius = radius;
}



// Definitions for PhysiksBase::Plane /////////////////////////////////////////

PhysiksBase::Plane::Plane(const Vector &normal, const Vector &origin) : ExtendedObject() {
    normal_vec = new Vector(normal.normalize());
    position.set(origin);
    unbounded = true;
}

PhysiksBase::Plane::Plane(const Vector &vec1, const Vector &vec2, const Vector &origin) : ExtendedObject() {
    basis_vec1 = new Vector(vec1);
    basis_vec2 = new Vector(vec2);
    normal_vec = new Vector(vec1.cross(vec2).normalize());
    position.set(origin);
    unbounded = false;
}

PhysiksBase::Plane::~Plane() {
    delete basis_vec1, basis_vec2, normal_vec;
}

void PhysiksBase::Plane::draw() const {
}

void PhysiksBase::Plane::bounded(bool flag) {
    unbounded = !flag;
}

bool PhysiksBase::Plane::bounded() const {
    return !unbounded;
}

bool PhysiksBase::Plane::contains(const Vector &vec) const {
// One would check for a solution to the system a*v_1 + b*v_2 = vec. What is the most efficient
// way to do so?
    return false;
}

Vector PhysiksBase::Plane::normal() const {
    return Vector(*normal_vec);
}

void PhysiksBase::Plane::set_normal(const Vector &norm) {
    normal_vec->set(norm.normalize());
}
/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja                                  *
 *   harish.kukreja@gmail.com                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/


using namespace std;
using namespace PhysiksBase;

int main(int argc, char *argv[]) {
    World *world1 = new World(3, 1.0E-5);
    Vector *pos1 = new Vector(3, 0.0);
    Vector *pos2 = new Vector(3, 0.0);
    Vector *vel1 = pos1;
    Vector *vel2 = pos2;
    pos2->set(2, 5.4234E6);
    PhysicalObject *obj1 = new PhysicalObject(45.54, 0.0, *pos1, *vel1);
    PhysicalObject *obj2 = new PhysicalObject(6.777E24, 0.0, *pos2, *vel2);
    world1->add_object(*obj1);
    world1->add_object(*obj2);
    world1->integrate(0.0, 5.0);
    
    delete pos1, pos2, world1, obj1, obj2;
    
    return 0;
}



/***************************************************************************
 *   Copyright (C) 2006 by Harish Kukreja                                  *
 *   harish.kukreja@gmail.com                                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//
// C++ Interface: world
//
// Description: Definitions for Physiks' 'World' class and its members
//
//
// Author: Harish Kukreja <harish.kukreja@gmail.com>, (C) 2006
//
// Copyright: See COPYING file that comes with this distribution
//
//


using PhysiksBase::World;
using PhysiksBase::PhysicalObject;
using namespace PhysiksBase::Computation;
using std::vector;
using std::cout;

// A whole bunch of overloaded World constructors -----------------------------
World::World(int dim, double tstep) :
        output(false), objects(0), force_fields(0), dimension(dim), timestep(tstep) {
}

World::World(int dim, double tstep, vector<PhysicalObject *> &init_obj) : output(false), objects(init_obj), force_fields(0), dimension(dim), timestep(tstep), time(0.0) {
}

World::World(int dim, double tstep, vector<PhysicalObject *> &init_obj, vector<ForceField *> &init_force_fields) : output(false), objects(init_obj), force_fields(init_force_fields), dimension(dim), timestep(tstep), time(0.0) {
}

//-----------------------------------------------------------------------------

World::~World() {}

void World::add_object(PhysicalObject &obj) {
    objects.push_back(&obj);
}

// Add a unary procedure with side-effects to the list of such procedures in the calling world.
void World::add_single_object_proc(void (*proc)(PhysicalObject *)) {
    single_object_procs.push_back(proc);
}

// Add a binary procedure with side-effects to the list of such procedures in the calling world.
void World::add_double_object_proc(void (*proc)(PhysicalObject *, PhysicalObject *)) {
    double_object_procs.push_back(proc);
}

void World::integrate(double start_time, double end_time) {
    this->time = start_time;
    while (this->time < end_time) {
        if (output) print_data();
        step();
    }
}

void World::step() {
    Computation::set_force_vectors_zero(objects);
    Computation::apply_force_fields(objects, force_fields);
    Computation::apply_acceleration_fields(objects, acceleration_fields);
    Computation::apply_velocity_fields(objects, velocity_fields);
    apply_user_defined_fields();
    apply_user_defined_procs();
    Computation::update_forces(objects);
    Computation::update_velocities(objects, timestep);
    Computation::update_positions(objects, timestep);

    this->time += timestep;
}

double World::current_time() const {
    return this->time;
}

PhysicalObject *World::object(unsigned int index) const {
    return objects[index];
}

unsigned int World::num_objects() const {
    return objects.size();
}

// Applies all user defined procedures passed to the world, to all active objects within the world.
void World::apply_user_defined_procs() {
    if (!single_object_procs.empty()) {
        for (uint i = 0; i < objects.size(); i++) (*single_object_procs[i])(objects[i]);
    }
    if (!double_object_procs.empty()) {
        for (uint i = 0; i < double_object_procs.size(); i++)
            for (uint j = 0; j < objects.size(); j++)
                for (uint k = j + 1; k < objects.size(); k++)
                    (*double_object_procs[i])(objects[j], objects[k]);
    }
}

void World::apply_user_defined_fields() {
    Vector temp(dimension, 0.0);
    if (!force_fields.empty())
        for (uint i = 0; i < force_fields.size(); i++)
            for (uint k = 0; k < objects.size(); k++) {
        temp = (*force_fields[i])(objects[k]->position);
        objects[k]->net_force += temp;
            }
// Incorporate electric and magnetic fields.
}

// Output needs to be clearer. Essential object distinguishers are omitted due to
// laziness.
void World::print_data() {
    for (uint i = 0; i < objects.size(); i++) {
        cout << time << ": " << objects[i]->net_force.core()->data[0] << ", " <<
                objects[i]->net_force.core()->data[1] << ", " <<
                objects[i]->net_force.core()->data[2] << "\n";
        cout << time << ": " << objects[i]->velocity.core()->data[0] << ", " <<
                objects[i]->velocity.core()->data[1] << ", " <<
                objects[i]->velocity.core()->data[2] << "\n";
        cout << time << ": " << objects[i]->position.core()->data[0] << ", " <<
                objects[i]->position.core()->data[1] << ", " <<
                objects[i]->position.core()->data[2] << "\n";
    }
}
